install.packages("dmetar")
install.packages("meta")
install.packages("metafor")
install.packages("corrplot")
install.packages("Hmisc")
library("dmetar")
library("meta")
library("metafor")
library("corrplot")
library("Hmisc")


GLPT1DMMETA_r= read.rm5(
  "GLP-1 for Type 1 Diabetes- metareg.rm5",
  sep = ",",
  quote = "\"",
  "GLP-1 for Type 1 Diabetes- metareg",
  numbers.in.labels = TRUE,
  debug = 0
)

c1o1 <- metacr(GLPT1DMMETA_r,comp.no=1, outcome.no=1)
year1 <- c(2016, 2016, 2016, 2019, 2018, 2016, 2019, 2020, 2018, 2015, 2020, 2012, 2013, 2020, 2020, 2016, 2016, 2016, 2016, 2016, 2016, 2020, 2014, 2021, 2022)
#Location:   MN = 1, NA = 2, Europe = 3, Asia = 4
loc1 <- c(1, 1, 1, 3, 2, 3, 3, 3, 2, 3, 1, 4, 4, 2, 3, 2, 2, 2, 1, 1, 1, 3, 2, 1, 3)
week1 <- c(26, 26, 26, 26, 52, 24, 52, 26, 24, 12, 26, 52, 52, 24, 26, 12, 12, 12, 26, 26, 26, 52, 52, 54, 12)
#Journal Article =1, Coference abstract =2 
pubtype1 <- c(1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)


#Regimen 1= lira 2= exena 3= albi
regimen1 <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 3, 2, 1, 1)
dose1 <- c(0.6, 1.2, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.2, 1.8, 0.9, NA, NA, NA, 0.6, 1.2, 1.8, 0.6, 1.2, 1.8, NA, NA, 1.8, 1.2)

#biases:   1= low risk, 2= medium, 3= high risk
rand1 <- c(1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 3, 3, 2, 1, 2, 2, 2, 1, 1, 1, 1, 3, 1, 2)
attr1 <- c(2, 2, 2, 3, 2, 2, 2, 1, 1, 3, 3, 2, 1, 2, 3, 3, 3, 3, 2, 2, 2, 1, 3, 3, 1)
blin1 <- c(1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1)
adhr1 <- c(2, 2, 2, 3, 2, 2, 2, 1, 1, 3, 3, 2, 1, 2, 3, 3, 3, 3, 2, 2, 2, 1, 3, 3, 1)
detec1 <- c(1, 1, 1, 1, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)


# cpep1 <- c(15.2, 15.2, 15.2, NA, NA, 21.0, 100, 9.1, NA, 0, NA, NA, 100, 41.7, 22, 0, 0, 0, 17.4, 17.4, 17.4, NA, NA, 100, 100)
# pump1 <- c(25.0, 26.0, 25.0, 17.9, NA, 0, NA, 100, 47, NA, 69, NA, NA, 73.4, 0, 80, 80, 80, 24.0, 28.3, 29.5, 0, NA, NA, NA)

age1 <- c(43.3, 42.8, 43, 50.4, 46.7, 48, 27, 46.5, 35.8, 37.8, 46, 48.5, 27.7, 36.1, 50.2, 47.5, 46, 46, 43.5, 43.7, 43.6, 22.5, 37.3, 29.5, 33.6) 
sex1 <- c(45, 47, 46, 80, NA, 65, NA, 32, 47, 66, 38, NA, 72, 32, 72.5, 45.7, 50, 34.4, 47.5, 48.2, 48.9, 56.5, 46, 65.3, 50)
dur1 <- c(20.9, 20.9, 21.1, 32.4, 22.3, 22.5, 0.1, 20.5, 20.3, 18.9, 18.6, NA, 0.1, 19.6, 21.1, 19, 20, 20, 21.3, 21.6, 21.6, 0.2, 20.5, 0.25, 3.4)
ins1 <- c(60.3, 60.3, 60.1, NA, NA, 59.5, 22.5, 51, 70.9, 59.6, 57.5, 39.2, NA, 49.5, 60.4, 49.5, 58.7, 47.1, 58.2, 58.5, 58.9, NA, 42.7, NA, 17.2)
bmi1 <- c(28.9, 28.9, 28.9, 29, 28.9, 30.1, 24, 29.5, 30.5, 23.4, 31.6, 22.8, 21.5, 29.3, 28.3, 27, 30.5, 28, 29.7, 29.6, 29.7, 22.3, 26.1, 24.1, 22.9)
a1c1 <- c(8.08, 8.08, 8.08, 8.2, 7.82, 8.7, 8.4, 8.15, 7.4, 8.75, 7.89, 7.65, 9.7, 7.6, 8.25, 7.57, 7.57, 7.57, 8.16, 8.16, 8.16, 7.29, 7, 7.25, 6.8)



factors1 <- data.frame(Publication_Year= year1,
                       Continent= loc1,
                       Study_Duration= week1,
                       GLP1_Subtype= regimen1,
                       GLP1_Dosage= dose1,
                       Publication_Type= pubtype1)
corm1 <- rcorr(as.matrix(factors1))
corm1$r[corm1$r < -2] <- NaN
corm1
corrplot(corm1$r, method = 'number')


factors1 <- data.frame(Randomization_Bias= rand1,
                       Attrition_Bias= attr1,
                       Performance_Bias= blin1,
                       Nonadherence_Bias= adhr1,
                       Detection_Bias= detec1)
corm1 <- rcorr(as.matrix(factors1))
corm1
corrplot(corm1$r,  method = 'number')


factors1 <- data.frame(Randomization_Bias= rand1,
                       Attrition_Bias= attr1,
                       Reporting_Bias= repo1)
corm1 <- rcorr(as.matrix(factors1))
corm1
corrplot(corm1$r,  method = 'number')



factors1 <- data.frame(Age= age1,
                       Male_Sex= sex1,
                       DM_duration= dur1,
                       Baseline_A1c= a1c1, 
                       Baseline_BMI= bmi1,
                       Daily_Insulin= ins1)
corm1 <- rcorr(as.matrix(factors1))
corm1
corrplot(corm1$r,  method = 'number')


factors1 <- data.frame(Male_Sex= sex1,
                       DM_duration= dur1,
                       Baseline_A1c= a1c1)
corm1 <- rcorr(as.matrix(factors1))
corrplot(corm1$r,  method = 'number')




factors1 <- data.frame(Publication_Year= year1,
                       Continent= loc1,
                       Study_Duration= week1,
                       GLP1_Subtype= regimen1,
                       GLP1_Dosage= dose1,
                       Publication_Type= pubtype1,
                       Randomization_Bias= rand1,
                       Attrition_Bias= attr1,
                       Male_Sex= sex1,
                       DM_Duration= dur1,
                       Baseline_A1c= a1c1)
corm1 <- rcorr(as.matrix(factors1))
corm1$r[corm1$r < -2] <- NaN
corm1
corrplot(corm1$r, method = 'number')



c1o1d= data.frame(TE= c(c1o1$TE), seTE= c(c1o1$seTE), 
                  Publication_Year= year1,
                  Continent= loc1,
                  Study_Duration= week1,
                  Publication_Type= pubtype1,
                  GLP1_Subtype= regimen1,
                  GLP1_Dosage= dose1,
                  Randomization_Bias= rand1,
                  Attrition_Bias= attr1,
                  Male_Sex= sex1,
                  DM_Duration= dur1,
                  Baseline_A1c= a1c1)

c1o1.mi <- multimodel.inference(TE= "TE", seTE = "seTE", data = c1o1d,
                                predictors = c("Publication_Year", "Continent", "Study_Duration",
                                               "GLP1_Dosage", "GLP1_Subtype", "Publication_Type",
                                               "Randomization_Bias",
                                               "Attrition_Bias",
                                               "Male_Sex", "DM_Duration", "Baseline_A1c"),
                                interaction=FALSE)
c1o1.mi
#
c1o1.mreg <- rma(yi = TE,
                 sei = seTE,
                 data = c1o1,
                 method = "ML",
                 mods = ~sex1,
                 test = "knha")
c1o1.mreg
c1o1.mreg1 <- rma(yi = TE,
                  sei = seTE,
                  data = c1o1,
                  method = "ML",
                  mods = ~dose1+sex1,
                  test = "knha")
c1o1.mreg1


anova(c1o1.mreg, c1o1.mreg1)



c1o2 <- metacr(GLPT1DMMETA_r,comp.no=1, outcome.no=2)
year12 <- c(2016, 2016, 2016, 2019, 2018, 2016, 2020, 2018, 2015, 2020, 2013, 2020, 2020, 2016, 2016, 2016, 2016, 2016, 2016, 2020, 2014, 2021, 2022)
#Location:   MN = 1, NA = 2, Europe = 3, Asia = 4
loc12 <- c(1, 1, 1, 3, 2, 3, 3, 2, 3, 1, 4, 2, 3, 2, 2, 2, 1, 1, 1, 3, 2, 1, 3)
week12 <- c(26, 26, 26, 26, 52, 24, 26, 24, 12, 26, 52, 24, 26, 12, 12, 12, 26, 26, 26, 52, 52, 54, 12)
pubtype12 <- c(1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)

#Regimen 1= lira 2= exena 3= albi
regimen12 <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 3, 2, 1, 1)
dose12 <- c(0.6, 1.2, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.2, 1.8, NA, NA, NA, 0.6, 1.2, 1.8, 0.6, 1.2, 1.8, NA, NA, 1.8, 1.2)


#biases:   1= low risk, 2= medium, 3= high risk
rand12 <- c(1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 3, 2, 1, 2, 2, 2, 1, 1, 1, 1, 3, 1, 2)
attr12 <- c(2, 2, 2, 3, 2, 2, 1, 1, 3, 3, 1, 2, 3, 3, 3, 3, 2, 2, 2, 1, 3, 3, 1)
detec12 <- c(1, 1, 1, 1, 2, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)


dur12 <- c(20.9, 20.9, 21.1, 32.4, 22.3, 22.5, 20.5, 20.3, 18.9, 18.6, 0.1, 19.6, 21.1, 19, 20, 20, 21.3, 21.6, 21.6, 0.2, 20.5, 0.25, 3.4)
a1c12 <- c(8.08, 8.08, 8.08, 8.2, 7.82, 8.7, 8.15, 7.4, 8.75, 7.89, 9.7, 7.6, 8.25, 7.57, 7.57, 7.57, 8.16, 8.16, 8.16, 7.29, 7, 7.25, 6.8)
sex12 <- c(45, 47, 46, 80, NA, 65, 32, 47, 66, 38, 72, 32, 72.5, 45.7, 50, 34.4, 47.5, 48.2, 48.9, 56.5, 46, 65.3, 50)


c1o2d= data.frame(TE= c(c1o2$TE), seTE= c(c1o2$seTE), 
                  Publication_Year= year12,
                  Continent= loc12,
                  Study_Duration= week12,
                  GLP1_Subtype= regimen12,
                  GLP1_Dosage= dose12,
                  Publication_Type= pubtype12,
                  Randomization_Bias= rand12,
                  Attrition_Bias= attr12,
                  Male_Sex= sex12,
                  DM_Duration= dur12,
                  Baseline_A1c= a1c12)

c1o2.mi <- multimodel.inference(TE= "TE", seTE = "seTE", data = c1o2d, method='REML', test='knha',
                                predictors = c("Publication_Year", "Continent", "Study_Duration",
                                               "GLP1_Dosage", "GLP1_Subtype", "Publication_Type",
                                               "Randomization_Bias",
                                               "Attrition_Bias", 
                                               "Male_Sex", "DM_Duration", "Baseline_A1c"),
                                interaction=FALSE)
c1o2.mi
